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1. INTRODUCTION 


This report summarizes the research activities accomplished 
during the 6 month period from 15 March through 14 September 1986. It 
includes a summary of work performed by personnel of Electro Magnetic 
Applications, Inc. (EMA) under their subcontract to this grant. 


2. WORK PERFORMED DURING THE REPORTING PERIOD 
2.1 IAS Research 

The work done during this period has centered on the completion of 
the first order parameterization scheme for the intracloud lightning 
discharge and its incorporation within the framework of the Storm 
Electrification Model (SEM). Previously, we established the criteria 
for the initiation, termination, and direction of propagation of the 
lightning channel path. This aspect of the parameterization scheme 
was incorporated into the framework of the SEM and tested in the con- 
text of the 19 July 1981 CCOPE simulation, with which we have been 
working. The aspect which remained to be dealt with was the redistri- 
bution of charge, created by the lightning, along the simulated chan- 
nel and, finally, into the space surrounding the channel. It is this 
charge which, when converted into ions, interacts with the other 
charges in the model and allows us to investigate the resultant 
effects of lightning on the subsequent electrical evolution of the 
model cloud. 

The charge redistribution problem was approached in the following 
manner. Using ideas from Kasemir (1964, 1980), we assume that the 
lightning channel is going to be electrically neutral over its entire 
length, indicating the deposition of equal amounts of positive and 
negative charge along the two portions of the simulated discharge. In 
the model, this charge deposition is accomplished by assuming that the 
charge density per unit length along the channel is proportional to 
the difference between the ambient potential, <(> (in Volts), at the 
grid point in question and the potential, $ 0 * the initiation point 
of the discharge. Thus 


Qo _ - k(^-<j» 0 ) 


( 1 ) 


where Q 0 is the linear charge density at a point along the channel, 
and k is a proportionality constant (farads/m). The negative sign 
assures that the charge deposited along that portion of the channel 
is basically opposite in sign to the ambient net charge in that grid 
volume, since the sign of the potential and the charge responsible for 
that potential are generally the same. 

The application of Eq. (1) results in a charge distribution along 
the channel similar to that shown in Fig. 1 between +1.2 and -1.5 km. 
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The two tails extending beyond these points will be explained 
subsequently. Since we are working with a discrete grid and the 
termination criterion forces the channel path to stop at a specific 
grid point (not at the exact place where this criterion is just met), 
charge neutrality is not guaranteed by this process. In order to 
overcome this difficulty, we arbitrarily extend the channel path four 
grid points beyond the designated termination point at each end. We 
first integrate over the channel path to find the total charge depos- 
ited on the sections of opposite polarity. We then compare these two 
and, at the end of the path segment having the greatest total charge 
magnitude, we cause an exponentially decreasing amount of charge to be 
deposited over the four additional grid points using the following 
expression 


Qj = Q n exp(-k{[j-n]A} 2 ) 


( 2 ) 


where j = n±l,...,n±4, Q n is the linear charge density (C/m) at the 
termination point of that segment [as determined by Eq. (1)], a is the 
grid spacing, and k is chosen such that Q n ±4 = Q n /1000. This new 
segment of charge density is integrated to obtain the total charge 
deposited along the segment of that polarity. We then calculate the 
absolute difference between the charge residing on the completed 
segment and the charge residing on the segment of opposite polarity. 
This difference is used as a basis for calculating a charge density 
decrease over the four grid points at the end of the remaining 
segment. 

As with the charge calculations on the initial channel, a 
trapezoidal integration is performed using Eq. (2) to represent the 
charge decrease over the four grid points. An iterative process 
(Newton-Raphson) is used on the integration scheme with the exponen- 
tial factor k being the unknown quantity. By varying k, the rate of 
decrease of charge density away from the initial endpoint of the seg- 
ment is controlled and the total additional charge distributed off the 
end of the segment can be set so that the total charge of that segment 
just balances the total charge on the other segment. By this proce- 
dure, the neutrality of the discharge channel is obtained over the 
discrete grid. 

Figure 1 shows the linear charge density calculated along the 
channel using Eq. (1) and a value of k = 2.5 x 10" 11 farads/m. The 
original channel is delineated by the 2 vertical bars on the horizon- 
tal axis. The effect of the charge balancing procedure is evident at 
the 2 tails. The total channel length is nearly 4.5 km. For this 
particular discharge (value of k), the total charge transferred is 
12 C. This may be a bit high for an intracloud discharge, however it 
is not an unreasonable value. 

To this point, we have the lightning discharge represented as a 
path on the discrete grid with a linear charge density at every grid 
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point along the path. In order to allow this discharge to interact 
with the other model parameters, we need to make some further modifi- 
cations. The discrete discharge path represents a line (ID) source in 
the 2D model. Because of the numerical methods employed in the model 
calculations, the sudden appearance of a line source of ions would 
represent a shock which the model could not accommodate. In order to 
make the channel ion production tractable for the model numerics, the 
charge must be distributed to some of the grid points adjacent to the 
channel path. A minimum of three grid points must be involved, the 
channel path itself and one grid point on each side. We have chosen 
to employ a total of 9 grid points (the path grid point plus four on 
each side) to represent the charge distribution. An expression iden- 
tical to Eq. (2) is used to specify the decay of the magnitude as a 
function of distance away from the channel, again with Q n ±4 = Q n /1000. 
For either a vertical or diagonal path orientation, the charge 
spreading is done in a horizontal direction. 

One additional consideration must be made. The charge density 
along the channel is expressed in terms of C/m. This must ultimately 
be converted into an ion density (ions/m 3 ). The linear charge den- 
sity, expressed by Q n , represents the charge density through a grid 
volume which surrounds the grid point and, therefore, can be converted 
to an equivalent volumetric charge density, assuming that the charge 
is distributed uniformly throughout the grid volume in question. In 
addition, when we undertake to distribute the line charge laterally 
away from the discharge path, we would be creating charge if we do not 
take steps to maintain charge conservation. However, charge conserva- 
tion can be accomplished during the spreading by adopting the 
following procedure. 

Figure 2 is a representation of the charge conservation concept. 

Q n represents the linear charge density at a particular grid point 
along the discharge path. Assuming that Qn represents the charge 
distribution in the grid volume AxAyAz surrounding the grid point, 
then the total charge in the grid volume is 


Qt = Qn Al 


(3) 


where a1 is the length of the discharge path segment for that grid 
volume. Because we are using a two-dimensional model, we must gener- 
ate a fictitious volume with which to work. Although, in a strict 
sense, the adoption of slab symmetry for the 2D model implies an infi 
nite extension to all model features in the y-direction, we choose 
here to adopt a concept of infinity which encompasses only the range 
of influence of the process being considered. Since we are employing 
a 9 grid point spreading procedure in the x-direction, we hypothesize 
a 9 grid point spreading in the y-direction also to comprise the 
influence volume. In Fig. 2, the rectangular volume at the center 
represents the total charge associated with the grid volume in 
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question. The two Gaussian type curves represent the redistribution 
of that charge over the 9x9 rectangular grid array surrounding the 
grid point where the total charge in both distributions is the same. 

Mathematically the procedure is developed as follows. The total 
charge under the two distributions is equated such that 


Al 4ay 4 ax _ k 2 _k V 2 

Qj = Q n Al * Q' J J ( e KlX e k 2 y dxdydz (4) 

o -4ay -4 ax 

which becomes 


4ay 4 ax 

Q t = Q n Al = 4 Q’a 1 J f e" kix2 e' k 2 y2 dxdy (5) 

o o 


which can be separated to yield 


4 £ y -kov2 4 £ x -k x2 

Qj = Q Al * 4Q a1 f e k 2 y dy f e k i x dx (6) 

o o 


Here, Q' is the magnitude of the charge density (C/m 3 ) at the grid 
point under consideration which will result in the same total charge, 
Q-f, when distributed under the exponential distribution rule. 

The two integrals are evaluated using 6 point Gauss-Legendre 
Quadrature with k x * k 2 , yielding a fall off to 1/1000 of the central 
point value at the 4th grid point [see discussion following Eq. (2)3. 
The result of this integration yields a relationship between Q' and Q n 
such that Q' * Q n /2.909 x 10 5 . Therefore, the linear charge density 
at a grid point may easily be converted to a volumetric charge density 
spread out laterally away from the discharge path in a manner identi- 
cal to Eq. (2), but with Q' replacing Q n . Figure 3 shows the results 
of the horizontal spreading of the charge distribution depicted in 
Fig. 2. Because of the exponential fall off, most of the charge is 
concentrated nearest to the main channel. The upper portion of the 
channel contains negative charge and the lower portion is positive. 

This process conserves the original total charge transferred by 
the discharge. Each calculated grid point charge density must be con- 
verted to an equivalent ion density. This is accomplished by dividing 
each charge magnitude by the electronic charge e = 1.6 x 10" 19 
C/electron assuming that all the ions produced are singly charged. 
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The ions resulting from the redistribution process are then added 
into the ambient ion fields as they exist at that point in the simula- 
tion. With the ion fields modified, a new total charge density is 
calculated, then Poisson's equation is employed to diagnose a new 
potential field from which the modified electric field components are 
calculated. This completes the discharge process and the time 
marching is resumed. The ions which were injected into the cloud 
along the discharge path interact with the charged hydrometeors in 
subsequent time steps and continue to modify the charge distribution 
and the resulting electric fields. 

The first experiment performed with the parameterization scheme 
used the value of k referred to above and resulted in a very strong 
influence on the subsequent electrification. Because of this, we 
decided to try a series of k values to see how differing amounts of 
charge transfer affected the electrical behavior. Using 10 values of 
k ranging from 2.5 x 10" 12 (weak charge transfer) to 2.5 x 10"H 
(strong charge transfer) we found a vastly different character to the 
electrical structure of the cloud. These results ranged from a minor 
perturbation for the weak cases to a complete domination of the charge 
structure and electrical evolution for the strongest cases. 

For demonstration purposes, some aspects of 4 of the cases have 
been chosen to illustrate these differences. In the following graphs, 
the numerical labels represent the following values of k [in Eq. (1)] 
and amounts of charge transferred: 

1 -> k = 2.5 x lO-n, 12 Coul 

2 -> k = 2.5 x 10-12, i.2 Coul 

5 -> k * 1.0 x 10-H, 4.8 Coul 

7 -> k = 1.5 x 10-n, 7.2 Coul. 

Figure 4 shows the variation of the maximum snow charge density (a), 
the maximum negative vertical (b) and maximum negative horizontal 
(c) electric field strengths for a period of 30 seconds following the 
time of the simulated lightning discharge for each of the four cases 
listed above. 

The majority of the snow in the simulated cloud is charged 
positively as a result of riming electrification. The lightning dis- 
charge occurs in the region of maximum snow charge density with both 
positive and negative ions being deposited amongst the positively 
charged snow. The negative lightning ions are formed in the region of 
maximally charged snow. Figure 4a shows that, in all four cases, the 
negative ions attach to the snow particles to reduce the maximum 
charge density. For the weakest case (2) the effect is very small and 
the electrical processes tend to return to their previous state. For 
the intermediate cases (5 and 7) there are sufficient negative ions 
for the attachment process to continue to erode the charge on the 
snow. For the strongest case (1) the trend is the same initially, but 
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then shows a large increase in the snow charge. Because the data used 
to make up this plot are the domain maximum values, this large 
increase represents the effect of the attachment of positive ions from 
the lower portion of the channel to snow particles lower in the cloud. 
While the original maximum charge region has been eroded by the 
negative ions, a new region of positive maximum has been created. 

Figure 4b shows the effect of the discharge ions on the maximum 
negative vertical electric field component. The weak and intermediate 
discharges all exhibit a tendency to reduce the field strength, with 
stronger discharges exhibiting the stronger tendency. This seems 
intuitively satisfying because we tend to think of lightning as acting 
to reduce the electrical stress in a cloud. For the strong case, how- 
ever, the results are somewhat unexpected. Here there is an initial 
i ncrease in the magnitude of the negative field component followed by 
a return to a near steady state at almost the initial value. Examina- 
tion of the structure of the vertical field component for cases 2, 5 
and 7 shows that the lightning causes an ever increasing perturbation 
on that structure; however, for case 1 an entirely different structure 
is established by the discharge. 

Figure 4c, which shows the maximum negative horizontal electric 
field component, further demonstrates the transition from perturbation 
to domination. Here, only the weak case shows a decrease in the field 
strength with a rapid return to the slow field buildup. The two 
intermediate cases both show an increase in the field strength imme- 
diately following the discharge; however, case 5 then shows a field 
reduction below the initial value, while case 7 continues to maintain 
field strengths greater than the initial value. The strong case shows 
a dramatic increase in this field component and maintains the 
increased magnitudes. 

While the parameterization scheme developed thus far is admittedly 
crude, it does seem to produce a lightning channel with realistic 
character. The parameter study indicates that the effects of a 
lightning discharge on the electrical evolution of a small thunder- 
storm depends quite heavily on the amount of charge transferred by the 
discharge. These effects range from being a minor perturbation on the 
system to dominating the electrical structure. For a vertical dis- 
charge, the most sensitive parameter to the lightning effects appears 
to be the horizontal electric field component. Because the actual 
nature of intracloud lightning is not well understood, a study of this 
type raises more questions than it answers. It is hoped that the stu- 
dies to be carried out during the remainder of this grant will lead to 
an improved understanding of the lightning process and thereby lead to 
an improved parameterization. 

2.2 EMA Subcontract - The Lightning Propagation Model 

Under subcontract to this grant, EMA has developed a Lightning 
Discharge Propagation Model (LDPM) for incorporation into the SEM. 
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The purpose of the LDPM is to account for the tortuous path that is 
usually followed by a lightning channel. One of the weaknesses of the 
SEM discharge parameterization scheme has been the propagation cri- 
teria which artificially constrain the channel to follow the side of a 
grid box or its diagonal depending on the ambient electric field dir- 
ection. The LDPM scheme takes into account the electric fields 
created by the leader tip and the previous leader segment. Also 
included are the effects of random free electrons in the vicinity 
of the leader tip as an influence in determining the direction of 
propagation of the next segment. 

The basic idea employed involves consideration of conditions in 
proximity to the leader tip. The leader tip contains a high concen- 
tration of charge. In the region around the tip, free electrons are 
randomly generated by the action of cosmic rays. These electrons are 
accelerated in the combined electric fields of the ambient atmosphere 
and the leader. If the acceleration is strong enough, such electrons 
can collide with and ionize neutral air molecules creating secondary 
free electrons which are subsequently influenced by the field. As the 
process proceeds, an avalanche of electrons will occur and the air 
will become highly conductive in the avalanche region. The LDPM 
algorithm determines which free electrons will result in the con- 
ducting path that will become the next leader step. By using random 
seeds, the LDPM produces lightning channels which propagate in the 
same general direction as the original scheme used in the SEM, but 
show deviations from the path more reminiscent of actual lightning. A 
report (EMA-86-R-55) covering the details of the LDPM is included as 
Appendix B. 

After completion of the LDPM development, a copy of the Fortran 
code was forwarded to us on magnetic tape. The code must be 
integrated into the current structure of the SEM. 


3. WORK TO BE ACCOMPLISHED 

In the time remaining under the grant, several tasks remain to be 
accomplished. The cloud-to-ground discharge parameterization scheme 
must be worked out (using the intracloud scheme as a basis) and incor- 
porated into the model. The LDPM must be worked into the model code 
and the termination criterion must be examined with consideration 
given to the electric field associated with the leader tip. A modifi- 
cation to the termination criterion is necessary because cloud-to- 
ground discharges are known to propagate through the subcloud region 
where the ambient field is often considerably below the threshold of 
150 kV/m. 

In addition, the details of the effects of the intracloud 
discharge on the electrification in the 19 July case will be examined 
because of the intriguing behavior exhibited during the variation of 
the charge transfer values in the parameter study. We think that this 
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approach is necessary because of the lack of current knowledge on the 
details of the lightning discharge. In short, what we have done thus 
far looks reasonable; however, we have no way to know that what we 
have done is right. 

Finally, the Wallops Island case and the 11 July 1978 TRIP case 
will be run with the lightning discharge active to see how the occur- 
rence of lightning, interacting with the other model variables, 
effects the subsequent electrification and what types of discharges 
will occur within the model. The 11 July TRIP case was chosen because 
it was a well documented, isolated storm and has been the subject of 
investigation by several other groups. 

The task remaining for EMA to accomplish under their subcontract 
involves the investigation of modifications to the lightning initia- 
tion criterion used in the SEM employing information they are 
acquiring on the effects of particle ensembles and their charges on 
the initiation of lightning. Because of the lack of field data on the 
locations and ambient conditions existent at the point of actual 
lightning, these studies will be speculative in nature. It is our 
hope that some new insights can be gained from this effort. 


4. CONFERENCE PAPERS AND PUBLICATIONS 

During the period covered by this report, one conference paper was 
submitted and subsequently accepted for oral presentation. The paper 
was entitled "A Lightning Parameterization Scheme in a Two-Dimensional, 
Time-Dependent Storm Electrification Model" authored by J. H. Helsdon, 
Jr., R. D. Farley, and G. Wu. It will be presented at the Conference 
on Cloud Physics sponsored by the American Meteorological Society to 
be held in Snowmass, CO, in September 1986. 

Also, a paper entitled "A Numerical Modeling Study of a Montana 
Thunderstorm: Part II. Model Results vs. Observations Involving 

Electrical Aspects" was completed and submitted for publication in the 
Journal of Geophysical Research . The paper was reviewed and accepted 
for publication after modification in light of the reviewer's com- 
ments. These modifications will be accomplished during the remainder 
of the grant period. Abstracts of the conference and journal papers 
are included in Appendix A. 

Finally, a Master's Thesis is being prepared by Mr. G. Wu. While 
Mr. Wu has not been supported directly by the NASA grant, the computa- 
tions that he has used in preparing his research have been supported 
by the grant. The topic of the thesis involves the lightning parame- 
terization scheme and its use in the SEM. Mr. Wu plans to finish his 
degree in time for December graduation. 
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Fig. 1: Linear charge density along the channel path using 

IqT~(T) and a value of k = 2.5 x 10" 11 f/m. The original channel 
is delineated by the vertical bars on the horizontal axis. 
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Fig, 2: Depiction of the 3 dimensional charge spreading procedure. 

The total charge in the rectangular volume and under the surface 
defined by the Gaussian curves are equal. 
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Fig. 3: Horizontal distribution of the linear charge in Fig. 1 after 

applying the procedure in Fig. 2. 
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Fig. 4: Domain maximum positive show charge density (a), maximum 

negative vertical electric field (b), and maximum negative horizontal 
electric field (c) for 30 sec following the discharge for four 
experiments. 
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APPENDIX A 


A LIGHTNING PARAMETERIZATION SCHEME IN A TWO-DIMENSIONAL, 
TIME-DEPENDENT STORM ELECTRIFICATION MODEL 

John H. Helsdon, Jr., Richard D. Farley, 
and Gang Wu 


ABSTRACT 


An important consideration in the modeling of thunderstorm 
electrification is the existence of lightning in nature. A model 
without lightning is only capable of simulating the development of a 
storm in its early stages (up to the time of first lightning). Since 
lightning accounts for significant charge rearrangement within a 
storm, it becomes necessary to include such processes in any modeling 
effort which hopes to shed light on the electrical evolution of such 
storms. In an attempt to make allowance for the lightning discharge 
in our Storm Electrification Model (SEM), we have developed a simpli- 
fied parameterization scheme to simulate the lightning channel and its 
attendant charge redistribution. 

This parameterization scheme has criteria for the initiation, 
propagation, termination and charge redistribution associated with 
lightning. In order to simplify the scheme as much as possible, the 
first three criteria are all dependent on the strength and direction 
of the electric field vector. The propagation of the channel is bi- 
directional and the charge distribution along the channel is calcu- 
lated so as to maintain overall charge neutrality for an intracloud 
discharge. The charge deposited along the channel is then spread in a 
conservative manner laterally away from the channel over several grid 
points in order to maintain computational stability. This distribu- 
tion of charges is then added to the small ion population and allowed 
to interact with the other charge carriers as the integration 
proceeds. 

A simulation of the 19 July 1981 CCOPE case has been run with the 
lightning discharge included. A parametric study has been undertaken 
to determine the effects of different amounts of charge transfer on 
the subsequent electrical evolution. The results of these simulations 
will be discussed. 


14 



A NUMERICAL MODELING STUDY OF A MONTANA THUNDERSTORM: 
PART II. MODEL RESULTS VS. OBSERVATIONS 
INVOLVING ELECTRICAL ASPECTS 

John H. Helsdon, Jr. and Richard D. Farley 


ABSTRACT 


In this investigation, we compare the results of the Storm 
Electrification Model (SEM) simulation of the 19 July CCOPE case study 
cloud against the actual observations with respect to the cloud's 
electrical characteristics, as deduced from the data of two aircraft. 

It is found that the SEM reproduces the basic charge and electric 
field structure of the cloud on a similar time scale to that observed. 
The comparability of the modeled and observed values of field strength 
and charge within the cloud is directly related to the proximity of 
the aircraft to the main active core of the cloud. The character of 
the electric field external to the cloud appears to be well modeled, 
at least at the time of observation. A feature which was not 
observed, but appears in the model, is a charge screening layer at the 
tip of the anvil . 

In addition, the model predicts electric field strengths 
sufficient to initiate a lightning discharge at approximately the same 
time in cloud evolution as that when an intracloud discharge was 
recorded. 
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A STOCHASTIC MODEL FOR THE TORTUOUS PROPAGATION 
OF A LIGHTNING STEPPED LEADER 

The results of Task 1 produced by EMA for the referenced contract 
"Completion of a Lightning Discharge Propagation Model (LDPM)" are discussed in 
this report. A FORTRAN code has been written for this model which can be integrated 
into the South Dakota School of Mines (SDSM) time-dependent Storm Electrification 
Model (SEM). The LDPM has been structured to enable its use in either a two or 
three-dimensional SEM with minor modifications. 

1. INTRODUCTION 

The LDPM is an approximation of the tortuous path of a lightning channel in 
a thunderstorm environment. Location and direction for the initial discharge point must 
be supplied as input data to the model. The discharge channel is then established as 
a series of connected stepped leader segments and propagates until certain 
termination conditions are met. All leader segments are assumed to be the same 
length and have equal charge densities. These parameters can be set to any 
values desired. While testing the LDPM, parameter values of 50 meters and 
1 millicoulomb/meter were used based on the average values published by Uman 1 . 

The primary parameter calculated in the model is the direction for a segment 
of the stepped leader. This computation is based on the effects of random free 
electrons in the vicinity of a leader segment tip which determine the direction of the 
next stepped leader segment. The details of the direction parameters for a stepped 
leader segment will be discussed in the next section. Once this direction is 
established, position vectors for a segment can be calculated. The model is then able 
to continue the computation of the positions and directions of following stepped leader 
segments. 


1 Uman. M.A.. " Lightning, ’ McGraw-Hill, New York, 1 969. 
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The remainder of this report discusses the calculations used in the LDPM, 
presents results obtained during the model tests, and provides a discussion on how 
the code is used. In addition, a listing of the code is provided. 

2. PHYSICAL CONSIDERATIONS 

Uman 1 discusses several theories which have been presented since the 
late 1 930’s to explain observed stepped leader characteristics. Although basic 
differences exist between these theories, several fundamental issues are common for 
the theories and are listed below. 

1 . A stepped leader has an inner conducting core wherein most of the 
current flows. 

2. The core is surrounded by a corona region in which the electric field 
exceeds the breakdown field of air (approximately 3 x 10 6 volts/meter 
for standard temperature and pressure at sea level). 

3. A mechanism exists which precedes the stepped leader to create an 
ionized channel in "virgin air" which provides continuation of the 
stepped leader propagation. 

4. As current flow increases into the ionized channel, a condition occurs 
whereby the conductivity increases until the channel becomes a new 
segment of the stepped leader and the end of the channel becomes the 
new tip of the leader. 

These theories agree that stepped leader segment directions are randomly 
distributed. Based on lightning channel tortuosity studies by Hill 2 * , the mean absolute 
direction change is nearly constant from flash-to-flash and on the order of 16 degrees. 
For purposes of the LDPM, it is assumed that the randomness of the stepped leader 
direction is related to the random creation of free electrons in air due to cosmic ray 
bombardment. 


2 Hill, R.D., "Analysis of the Irregular Paths of Lightning Channels," J. Geoohvs. Res .. Vol. 73, 1 968 

(pp. 1897-1906). 
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Issue 3 is of primary concern for this model. The stepped leader 
propagation theories often vary when describing this mechanism. The process 
occurring which results in the formation of the ionization channel ahead of the 
stepped leader tip has been more recently described by Rudolph 3 from studies 
involving the interaction of lightning with aircraft in flight. 

The leader tip is a region which contains a large charge density and high 
electric field levels. Free electrons (generated by cosmic rays) in the air are 
accelerated in this high field until they collide with a neutral atom or molecule. If the 
electron's kinetic energy is large enough at the time of the collision, the neutral particle 
can have an electron separated from it, producing a second free electron and a 
positive ion. The free electrons are then again accelerated by the field, possibly 
suffering more collisions and producing more free electrons and ions. If the rate of 
production of free electrons is larger than the rate of loss (by recombination, and 
attachment to form negative ions), an electron "avalanche" occurs, in which sufficient 
numbers of electrons and ions are produced and substantially alter the electrical 
conductivity of the air. 

Figures 1 and 2 illustrate how the effects of avalanching result in the 
generation of a new leader step from an existing leader tip. Figures 1 and 2 also 
illustrate the ionization path mechanisms for a negatively charged and for a positively 
charged leader tip, respectively. 

In Figure 1 , a free electron (at time ti) is accelerated by the electric field 
towards the leader tip. As avalanching occurs (time t 2 ), a net positive charge is formed 
in the region of the initial free electron position. This results from positive ions which 
have a much lower mobility or drift velocity than electrons. The net effect is 
neutralization of the electric field at the leader tip (time t 3 ). The conductivity in the 

region approaches the conductivity in the inner core of the stepped leader. Thus, the 
region in which positive ions were first produced becomes the new leader tip. 

In Figure 2, a free electron (at time t-|) is accelerated by the electric field 
away from the leader tip. As avalanching occurs (time t 2 ), a net positive charge is 

3 Rudolph, T„ and R.A. Perala, “Linear and Non-Linear Interpretation of the Direct Strike Lightning 
Response of the NASA FI 06B Thunderstorm Research Aircraft," Electro Magnetic Applications - Denver, 

EMA-83-R-21, March 21, 1983. 


3 



4 


Figure 1. Time Scenario of Free Electron Avalanching for 
Positively Charged Leader Tip 



9 

o 



5 


Figure 2. Time Scenario of Free Electron Avalanching for 
Negatively Charged Leader Tip 



formed in the region of the initial position of the free electron. Moreover, as the 
electrons move away from this region, a conducting channel is generated between the 
positive ions and electrons until the electric field between the two opposite charge 
regions is neutralized. However, the electric field is enhanced between the leader tip 
and the positive charge region and induces additional forces on the negative charge 
region around the tip. In addition, the breakdown region is extended toward the 
positive charge region. This effect is illustrated in Figure 2 at time t 3 as a new 

conductive region extends the stepped leader. The leader tip occurs at the negative 
charge region created by the electrons which were first accelerated away from the 
original free electron position. 

The mechanism for extending the stepped leader is dependent on the 
location of free electrons which create a conductive channel to the leader tip. Thus, it 
is necessary to consider all free electrons and determine which one will initiate the first 
complete ionized path to the leader tip. Moreover, the algorithm used must determine 
which electron will initiate this complete ionized path in the minimum time 
(compared to other electrons). 

A rigorous mathematical treatment of the mechanisms just discussed would 
require involved calculations for the avalanche process and resultant conductivity in 
the stepped leader tip vicinity. These calculations would also involve significant 
computer time when used in the SEM. A simplified algorithm which has similar 
dependencies on the position of random electrons and the electric field around the tip 
of a stepped leader segment was developed. 

Since the drift velocities of the positive ions are significantly lower than those 
of electrons, the first simplification does not account for the mobility of resultant positive 
ions. The equation of motion for electrons generated by the avalanche is then given 
by: 


dV e 
cf It 


qe_ 

m e 


(E T + V e x B ) - 



0 ) 


where: q e 

m e 


electron charge, 
electron mass, 
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Ej = total electric field, 

V e = electron velocity, 

B = magnetic field, and 
v = an approximate collision frequency. 

To further simplify this equation, it is assumed that the creation of electrons 
due to avalanching has progressed to a level where a collision dominated, plasma 
region exists. In turn, the electron velocity is constant and equal to the drift velocity, 
and: 



= 0 . 


(2) 


Thus, 

^-E-»VO, or V d = • MeE (3) 

where: (Iq = electron mobility, and 

V d = drift velocity. 

It is assumed that the vector r e between the leader tip and the electron 
position can be approximated by: 

r e = V d t, (4) 

or the time necessary to create a channel can be approximated: 


t = — oc — ^ . (5) 

IVdl Ie t | 

It is assumed in the model that the ionized channel induced by the 
avalanching process (resulting from randomly occurring electrons) will initiate the next 
stepped leader segment in the same direction for either a positively or negatively 
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charged stepped leader segment. Thus, the model is independent of the leader tip 
charge. Based on these assumptions, the new stepped leader segment direction is 
determined by identifying the free electron in the vicinity of the previous stepped 
leader segment tip which will create an ionization path in a minimum time (compared 
to other free electrons). 

The model requires the following computations which will be discussed in 
detail in the next sections. 

a. Random generation of electrons in the vicinity of each stepped leader 
tip. 

b. Total electric field strength and direction at each electron location. 

c. The direction vector of the new stepped leader. 

3. LIGHTNING DISCHARGE PROPAGATION MODEL 

CALCULATIONS 

3.1 Description 

The LDPM code was written to predict the tortuous path of a lightning 
channel based on the physical considerations in Section 2. The routine is 
programmed to calculate the parameters for each stepped leader segment. 

3.1.1 Random Electron Generation 

A Monte Carlo technique was used to generate the random locations of free 
electrons within a sphere around each leader segment tip. The ambient electron 
density is assumed to be .2 electrons/m 3 . The size of the sphere is dependent upon 
the number of free electrons. Each model run in this report considered 500 random 
electrons around each tip. Thus, the sphere centered at the leader tip was required to 
have a radius of 8.419 meters. 
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A seed integer is also required for the Monte Carlo routine. This seed 
controls the random placement of electrons. Different seeds will produce different sets 
of electron locations. 

3.1.2 Total Electric Field Calculation 

The total electric field strength and direction was computed for each electron 
location generated by the Monte Carlo technique. The electric field at any electron 
location should be the resultant of the ambient field and the fields due to the previous 
leader segments. The model code was tested to incorporate the effects of the previous 
leader segments and the ambient field in the electric field calculation. It was observed 
that this calculation only required the ambient electric field and the electric field 
components due to the two leader segments closest to the electron. Electric fields 
caused by segments further away did not affect the modeled lightning path 
significantly. The total resultant electric field at a location (x,y,z) in the vicinity of the ith 
stepped leader segment is stated below: 

E total ( x .y.z) = E ambient ( x .y.z) + E i segment (x.y.z) + E i -1 segment (x,y,z) . (6) 

The general equation for the electric field at any point in space due to a line 
charge is given below and is illustrated by Figure 3. 

The electric field at any point in space due to the line charge element Xdl is: 


dE 


Xdl (r - r') 
4tce 0 | r- ?|3 


; X = charge/unit length . 


Thus, the electric field at point p due to the total line charge is: 


E 


X 

47CE 0 



r - (R + lu) 
r-(R + lu) 1 3 


( 7 ) 


( 8 ) 
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Figure 3. Stepped Leader Segment Line Charge and 
Its Parameters in World Coordinates 


The general world coordinate solution for the electric field due to the line charge 
becomes: 


X 

A(L * B) + u(A2-BL) 


A 2 u - AB 


4 to 0 (A2 - B2) j 

(A2- 2BL + L2)1/2 


|Al 



where: A = r- R , and 

B = A • u. 
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Figure 4 depicts the constant electric field contours due to a line charge. 
The arrows show the direction of the electric field 


3.1.3 Direction Vector of the New Stepped Leader Segment 

Only one of the electrons generated in the vicinity of the leader tip will 
properly satisfy the model requirements and determine the direction of the new 
stepped leader segment. For simplification, each stepped leader segment is assumed 
to be positively charged (as discussed in Section 2) for the LDPM. The new stepped 
leader segment is assumed to be initiated by the ionized path formed in the least time. 
The time to generate the ionized path is proportional to the ratio formed from the 
distance of the leader tip to the electron (r e ) and the resultant electric field strength (Ej) 

at the electron location. This relationship was determined in Section 2 and is stated in 
Equation 5. Thus, the new stepped leader segment direction is chosen to be the 
electric field direction at the initial location of the free electron which initiates an 
ionized path in the minimum time (compared to other electrons). 


If the angle 0 e , formed by the Ej and r e vectors at the initial electron location 
is within a certain angle range, that electron will create an ionized channel to the 
leader tip along a "relatively" straight path. The time it takes to establish the ionization 
channel can be approximated using Equation (5). As 0 e increases beyond this certain 

angle range, the actual path of the ionization channel along electric field lines will 
deviate from a straight line. Therefore, the actual time it takes to establish the ionized 
channel to the leader tip may be significantly greater than the time approximated by 
Equation (5). In fact, if 0 e is large enough, the ionized path may not attach to the 
leader segment tip, but to some point behind the tip instead. The cosine of 0 e formed 
by Ej and r e at any electron location is computed as follows: 


cos 0 e 


4 * E t 
l?.l I Erl 


( 10 ) 


A stepped leader segment is defined by the computer code as a finite length 
with its tip being a single point. In nature, the leader tip is not a single point but is 
actually defined as a volume. In terms of the LDPM, this volume exists around the 
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Figure 4. Constant Electric Field Amplitude Contours for Line Charge 
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leader tip point and has an electric field strength greater than 3 x 10 6 V/m everywhere 
inside. 


Figure 5 depicts the above concepts. The angle range parameter is 
determined by the radius of the corona region. This radius is on the order of 2 meters 
and is assumed to be the same magnitude for all leader segments in the LDPM. The 
angle range parameter was adjusted to an absolute value of 1 2.5 degrees. This value 
appeared to allow the lightning channel path to be reasonably tortuous. The angle 
range parameter can be changed in the LDPM as desired. However, 1 2.5 degrees is 
recommended, since smaller angle parameter values limit the tortuousity, while larger 
angle parameter values result in erratic tortuousity which exceeds the findings of Hill 2 . 

Each of the 500 electrons generated by the Monte Carlo routine is tested in 
the code to determine which one is used to establish the direction of the new stepped 
leader segment. The one electron chosen must meet the following conditions: 

1 . The initial free electron which initiates an avalanche must originate in a 
location where the electric field strength is less than 3 x 10 6 V/m. (If the 
field was larger, the electron was assumed to be inside the corona 
region of the leader segment and was ngi considered). 

2. The angle 0 e formed by Ex and r e for an electron location must be less 

than the angle range parameter (absolute value of 1 2.5 degrees) 

(If 0 e was greater, the electron was uol considered). 

3. Of the remaining electrons, the electron which creates an ionization 
channel to the leader tip in the least time, t (as computed in 
Equation (5)) was identified and used to calculate the new direction. 

The new stepped leader segment direction is then taken to be the direction of the total 
electric field at this "identified" electron location. 

It is important to note that, if the Monte Carlo sphere is too small, most or all 
of the random electrons will lie in a region where the electric field strength is above 
3 x 1 0 6 V/m. In turn, the model may not function properly. It was observed that a 
sphere containing 500 electrons and having a radius of 8.419 meters worked best 
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Electron "C" will not be drawn to the leader tip at an and equation [5] does not apply. 

However, 9 e >12.5* and the electron will not be considered. 

Figure 5. Paths of Free Electrons in Vicinity of Stepped Leader Tip 
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during the LDPM development. This sphere size contained enough electrons in a 
region where the electric field was below breakdown and allowed the LDPM to 
function properly and efficiently. 

3.2 Model Termination Conditions 

The LDPM will terminate if any of the following conditions are met: 

1. The stepped leader enters a defined region (grid cell) with an ambient 
electric field less than 1.5 x 10 5 V/m. This termination criterion was 
based on work done by Griffiths 4 involving the propagation of positive 
streamers in the laboratory, 

2. The stepped leader propagates outside the designated problem space, 

3. (Optional) The stepped leader enters a defined charged region 
(Specified grid cell). 

3.3 MODEL INPUT AND OUTPUT REQUIREMENTS 

i 

3.3.1 Input 

The problem space is depicted as a rectangular grid either in two or three 
dimensions. The code is written to accommodate two or three-dimensional data. This 
version of the LDPM is adjusted for a two-dimensional SEM. Each grid cell has side 
lengths equal to 200 meters. This routine requires the number of grid cells along the 
x-axis (horizontal), along the z-axis (vertical), and the initial discharge coordinate (x,z). 
Ambient electric field components are required for each grid cell. The routine also 
requires a seed integer for the Monte Carlo scheme. These variables are: 

1. NCELX Number of cells along the x-axis 

NCELZ Number of cells along the z-axis, 


4 Griffiths, R.F. and C.T. Phelps, "The Effects of Air Pressure and Water Vapor Content on the Propagation 
of Positive Corona Streamers and Their Implications to Lightning Initiation," Quart. J. Rov. Meteor. Soc .. Vol 

102, 1976 (p. 419). 
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2 . 


ICX, ICZ 


Initial discharge point (x,z), 


3. FEX(K), FEZ(K) One-dimensional array of x,z ambient field 

components, and 

4. ISEED Monte Carlo seed integer; must be less 

than 32,000. 

For the NCAR Electricity Model: 

1. NCELX = 97; NCELZ = 97 

2. ICX = 86; ICZ = 31 

3. FEX(K), K = 1 , 9409; FEZ(K), K = 1 , 9409 

3.3.2 Output 

The routine generates the (x,z) position in world coordinates for each 
stepped leader segment tip in the propagation path. The x and z arguments are 
written in the one-dimensional RX and RZ arrays. Coordinate values x and z are given 
in meters from the origin (x = 0, z = 0). These arguments can be passed back into the 
main program and are as follows: 

1. RX(NL), RZ(NL) One-dimensional arrays of the x and z values in 

meters 

NL corresponds to the leader number (not the 
grid cell number), and 

2. NFL Integer; total number of NL leader segments. 

3.4 MODEL EXAMPLES 

Several lightning discharge propagation paths were generated and are 
shown in Figures 6 through 8. Figure 6 shows the lightning paths modeled for the 
NCAR electricity data at 58.5 minutes. The solid and dashed lines are the models 
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Figure 7. EMA Electricity Model Data - Initial Discharge 
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(Comparison of Monte Carlo Seeds) 
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produced from the LDPM and the algorithm described in Helsdon 3 * 5 , respectively. The 
four plots in Figure 6 demonstrate the effects of varying ISEED (Monte Carlo Seed for 
Generating Random Electron Locations). 

Figure 7 shows the effects of ISEED for ambient electric fields generated 
from an in-house routine. This routine creates electric fields from a set of discrete 
charge centers. Figure 8 shows the effects of varying the charge magnitudes at 
specified cells. These ambient fields were computed from the in-house routine. 

4. MODEL VERSIONS AND DATA 

The listings for three versions of LDPM are included in Appendix A of this 
report. In addition, the code used to generate ambient electric fields is included. The 
listings are as follows: 

1 . LDPM - Subroutine 

Subroutine with common statements; requires the number of x cells 
(NCELX), z cells (NCELZ), the initial discharge coordinate, ambient 
electric field components at each cell, and ISEED. Passes leader tip 
locations RX(I), RZ(I) and NFL, the number of total leaders back into 
main program. 

2. LDPM.FOR 

Routine reads in ambient electric fields; requires NCELX, NCELZ, 
initial discharge coordinate, and ISEED. Calculates leader tip 
locations. 

3. LDPM - Interactive (not required by contract; used as test program) 
Interactive routine; generates ambient electric field components over 
grid from specified cells and their charge magnitudes. Requires 
NCELX, NCELZ, the initial discharge coordinate, and ISEED. Routine 
also requires coordinates of specified charge cells. Calculates leader 

tip locations. 

5 Helsdon, J.R, ft.b. Farley and G. Wu, "A Lightning Parameterization Scheme in a Two- 
Dimensional, Time-Dependent Storm Electrification Model," presented at the 1986 Conference 
on Cloud Physics, held Sept. 22-26, 1986 in Snowmass, Colorado. 
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4. EFS.DT 

NCAR ambient electric field data (58.5 minutes, unformatted) 
FEX(K), K = 1 , 9409; FEZ(K), K = 1 , 9409. 


21 


ADDITIONAL REFERENCES 


1. Krider, E.P., C.D. Weidman, and R.C. Noggle, "The Electric Fields Produced by 
Lightning Stepped Leaders", J. Geophv. Res .. 82, 1981 (pp.951-960). 

2. Fleagle, R.G., J.A. Busingeil, " An Introduction to Atmospheric Physics ". 
Academic Press, New York, 1963. 

3. Golde, R.H., " Physics of Lightning - Volume 1 Academic Press, New York, 
1972. 

4. Hill, R.D., "Analysis of Irregular Paths of Lightning Channels," J. Geophvs. Res .. 
Vol. 73, 1968 (pp. 1897- 1906). 


22 


APPENDIX A 
LDPM CODE VERSIONS 


A-1 



LDPM -SUBROUTINE 
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SUBROUTINE LDPM 

DEVELOPED BY E.M.A. - AUGUST 1986 

SUB-PROGRAM DETERMINES DIRECTION OF NEXT LEADER SEGMENT DUE TO 
ELECTRON WHICH WILL FORM AN IONIZED PATH TO 
THE STEPPED LEADER SEGMENT TIP IN THE LEAST TIME 
E FIELD TOTAL * E FIELD LEADER + E FIELD PREVIOUS LEADER 
+ AMBIENT E FIELD 


NCEL= * CELLS IN X AND IN Z 

ICXfICZ* RXrRZ COORD OF INITIAL CHARGE REGION 

NFL* TOTAL NUMBER OF STEPPED LEADER SEGMENTS 

C0MM0N/MAIN/RX<O:500) fRZ'OJSOO) fFEX<9409> fFEZ(9409) 

COMMON/V AR/NCELX f NCELZ fICXf ICZfNFL 
CHARACTER*^ 0FILE1fQFILE2 

ISEED=RANDOM FUNCTION INTEGER 
ISEED=30?S7 

QFILE1= / LDPM*DAT' ! BINARY PLOTFILE OF LIGHTNING PATH 
OFILE2='LDPM*PRINT' ! PRINT VALUES OF LIGHTNING MODEL 

INPUT VARIABLES 

CL*SEGMENT LENGTH ELAMBDA=LEADER CHARGE DENSITY 
BDF=BREAK DOWN FIELD EK=1/4PIEpso NL=LEADER * 

DZ=Z CELL HEIGHT DX=X CELL WIDTH 

NE=4 FREE ELECTRONS ED=ELECTRON DENSITY U2/»**3> 

Q=MAX ALLOWABLE ANGLE BTW EFIELD l DIRECTION VECTORS 

DATA CL f EL AMBDA f BDF/50 ♦ f ♦ 001 f 1 , 5E5/ 

DATA PIfEKfNL/3* 1415927f9*E9f0/ 

DATA DZfDXfNEfEDf Q/ 200 * f 200 ,» 500 1 ♦ 2» 12 . 5/ 

DATA EXM1fEYM1fEZM1/3*0»/ 

TWQPI*2**PI 

QP*Q*PI/180* 

CONST =EK*ELAMBDA 
VOL =FLO AT ( NE ) /ED 
RADIUS* ( * 75*V0L/PI ) ** < 1 * /3 . ) 

KTOT*NCELX*NCELZ 

XTOT=NCELX*DX 

INITIALIZE STEPPED LEADER PROPAGATION 
INPUT DATA: Rl, AMBIENT E FIELBf U1 
COORDINATES OF INITIAL POINT ON GRID PROBLEM SPACE 
RX<0)-DX* (FLOAT < ICX)-*5) 

RY=0. ! FOR 2D CASE? RY(0)*DY*FL0AT<ICY)-*5> r FOR 3D CASE 

RZ ( 0 ) =DZ* ( FLOAT ( ICZ) -*5 ) 


IC*ICX+NCELX*<ICZ-1) 

FOR 2D CASE 

ETQT=SGRT ( FEX ( I C ) *FEX C I C )+FEZ ( IC ) *FEZ < IC ) ) 

FOR 3D CASE 

ETOT =SQRT ( FEX ( IC) tFEX ( IC ) +FEY ( IC ) *FEY ( IC ) +FEZ ( I C ) tFEZ ( IC ) ) 
DEFINE UNIT VECTOR 
UX=FEX(IC)/ETOT 

UY*0# ! FOR 2D CASE 5 UY=FEY<IO/ETOTf FOR 3D CASE 

UZ=FEZ( IO/ETQT 

GPEN< 1 fFILE=0FILE1fF0RM*'UNF0RHATTED' fSTATUS='NEW') 

OPEN ( 2 fFILE=QF I LE2 f STATUS* / NEW / ) 
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C START LOOP 
2000 NL=NL+1 
C 

C DETERMINE CELL POSITION OF NEW LEADER SEGMENT TIP ON OVERLAY 

C USE APPROPRIATE AMBIENT E FIELD MAGNITUDE AND DIRECTION 

C K IS THE CELL NUMBER OF THE LEADER SEGMENT TIP LOCATION 
KX=ANINT ( (RX(NL-l)/DX)+.5) 

KZ=ANINT( (RZ(NL-1 >/DZ)+.5> 

K=KX +NCELX* ( KZ-1 ) ! NEED KY FOR 3D CASE 
C FOR 2D CASE 

EFAMB=SQRT(FEX(K) *FEX ( K HFEZ ( K ) *FEZ ( K) ) 

C FOR 3D CASE 

C EFAMB=SQRT < FEX < K> *FEX ( K HFE Y < K > *FEY < K> +FEZ ( K ) *FEZ ( K ) ) 

C 

C TERMINATE IF 

C LEADER SEGMENT HAS PROPAGATED OUTSIDE THE GRID 
C LEADER SEGMENT ENTERING CELL KITH AN AMBIENT FIELD LT 1.5M V/M 
IFIK.LT.l .OR.K.GT .KTOT .OR. 

+ RX ( NL-1 ) . LT . 0 .OR . RX ( NL-1 ) . GT .XTOT ) GOTO 3000 

C 

URITE(2»*> 'LEADER SEGMENT NUMBER' »NL 

URITE(2» *) 'CELL NUMBER ',K»' <X»Z> COORD'.KX.KZ 

URITE(2.*)'T0T E FIELD '.ETOT. ' AMB E FIELD ' »EFAMB 

URITE(2.*)'AMB EX '.FEXOO.' AMB EZ '.FEZOO 

URITE(2.*)'X COORD ' .RX(NL-l) .' Z COORD '»RZ(NL-1> 

URITE<2»*>' 

UR I TE ( 1 ) RX ( NL-1 > » RZ ( NL-1 ) 

C 

IF ( EF AMB » LT . BDF ) GOTO 3000 
C 

C SHORTEST TIME OF ELECTRON TO FORM IONIZED CHANNEL 
C THIS TIME DETERMINES DIRECTION OF NEU LEADER SEGMENT 
TMIN=1.E20 
C 

C USE MONTE CARLO METHOD TO GENERATE RANDOM FREE ELECTRONS 

C TO CALCULATE ELECTRON DISTRIBUTION NEAR LEADER SEGMENT TIP 

RXP=RX(NL-1)+UX*CL 

RYP=RY+UY*CL ! USE RY(NL-l) FOR 3D CASE 
RZP=RZ<NL-mUZ*CL 
C 

DO 100 1=1. NE 

IN UORLD COORDINATES 

CALL RANDOM (I SEED. RAND) 

RE=RAND*RADIUS 

C RE=SQRT ( RANDtRADIUSftRADIUS ) ! FOR ACTUAL RANDOM DISTRIBUTION 

CALL RANDOM ( ISEED . RAND ) 

TH=RAND*PI 

C TH=AC0S((RAND*2)-1) ! FOR ACTUAL RANDOM DISTRIBUTION 

CALL RANDOM ( ISEED » RAND ) 

PH=RAND*TUOPI 
RMX=RE*SIN ( THXCOS < PH ) 

RMY=RE*SIN(TH)*SIN<PH) 

RMZ=RE*COS(TH) 

RM=SQRT( RMXTRMX+RM Y*RMY+RMZ*RMZ ) 

X=RMX+RXP 

Y=RMY+RYP 

Z=RMZ+RZP 

C 

C COMPUTE E FIELD COMPONENTS FOR Ith ELECTRON 
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AX=0C-RX(NL-1) 

AY=Y-RY ! USE RY(NL-l) FOR 3D CASE 
AZ=Z-RZ(NL-1> 

A-SQRT ( AX* AX+ AY*AY+AZ*AZ ) 

IF<A*LT.1*E-4)GGT0 100 

B=AX*UX+AY*UY+AZ*UZ 

A5-A*A 

BS=B*B 

IF'AS-BS.LT.1*E-6)GGT0 100 
P=CQNST / ( AS-BS > 

D-CL*CL-2 . *B*CL+AS 
IF(D*LE»0* )GOTO 100 
D=SGRT(D) 

C 

C E FIELD COMPONENTS DUE TO LEADER SEGMENT AND AMBIENT FIELD 
EX=P*< ( ( AX*<CL-B)+UX*(AS-B*CL) >/D)- 
+ ( <UX*AS-AX*B)/A) )+FEX(K> 

EY=P*( < < AY*(CL-BHUY*(AS-B*CL) )/D)- 
+ <<UY*AS-AY*B)/A)> l+FEYCK) FOR 3D CASE 

EZ=P*< ( ( AZ*<CL-B)+UZ*( AS-B*CL) )/D>- 
+ ( ( UZ*AS- AZ*B )/A) > +FEZ ( K ) 

C 

IF(NL*LE*1)G0TQ 1000 
AXMl=X-RX<NL-2> 

AYM1=Y-RYM1 ! FOR 2D CASE? RYMl=RY(NL-2>* FOR 3D CASE 
AZMl=Z-RZ(NL-2) 

AM1=SQRT(AXM1*AXM1+AYH1*AYM1+AZM1*AZM1> 

IF(AM1 *LT * 1 fE-6)G0T0 100 

BM 1 =AXM 1 *UXM1 +AYM1*UYM1+AZM 1 *UZM1 

ASM1=AM1*AM1 

BSM1~BM1*BM1 

IF( ASM1-BSM1 *LT * 1 *E-6)G0T0 100 
PM 1 =CQNST / ( ASM1 -BSM 1 ) 

DM1 =CL*CL-2 * *BM1*CL+ASM1 
IF(DM1 »LE»0* )G0T0 100 
BM1=SGRT(DM1> 

C 

C E FIELD COMPONENTS DUE TO PREVIOUS LEADER SEGMENT 

EXM1=PM1*<<<AXM1*<CL-BM1)+UXM1*<ASM1-BN1*CL)>/DM1>- 
+ < ( UXM1*ASH1-AXM1*BM1 ) /AMI ) ) 

EYM1=PM1*< ( (AYM1*<CL-BN1)+UYM1*(ASM1-BM1*CL) )/DMD- 
+ ( <UYM1*ASM1-AYM1*BM1)/AM1 ) ) 

EZM1=PM1*( ( < AZM1*(CL-BM1 )+UZMl*(ASMl-BMl*CL) )/DMl)- 
+ ( (UZM1*ASH1-AZM1*BM1)/AM1 ) ) 

C 

C TOTAL E FIELD COMPONENTS FOR Ith ELECTRON 
1000 EX=EX+EXM1 
EY=EY+EYM1 
EZ=EZ+EZM1 

ET=SORT ( EX*EX+EY*EY+EZ*EZ H * 01 
DL= ( RMX*EX+RMY*EY+RMZ*EZ ) / ( RM*ET ) 

C 

C CHOOSE APPROPRIATE ELECTRON TO FORM IONIZED CHANNEL 
C TO THE LEADER SEGMENT TIP IN SHORTEST TIME 
IF(ET *LT* (3.0E6) *AND*(RM/ET) •LT.TMIN* 

+ AND * ACOS ( ABS < DL ) ) * LE * QP ) THEN 
THIN=RM/ET 
ETOT=ET 
DIRX=EX 
DIRY=EY 
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DIRZ-EZ 
END IF 

100 CONTINUE 
C 

C CALCULATE DIRECTION FOR NEXT LEADER SE6MENT 
DRX=CL*UX 
DRY=CL*UY 
DRZ=CL*UZ 
RYM1=RY 

RY=RY+DRY ! FOR 2D CASE* R Y ( NL ) =RY ( NL-1 ) +DRY ? FOR 3D CASE 
RX ( NL ) =RX ( NL-1 ) +DRX 
RZ ( NL ) =RZ ( NL-1 ) +DRZ 
C 

C NEW UNIT VECTOR 

C UNIT VECTOR IS PROJECTED ONTO XfZ PLANE FOR 2D SPACE 
C FOR 3D CASE 

C DIR=SQRT<DIRX*DIRX+DIRY*DIRY+DIRZ*DIRZ) 

C FOR 2D CASE 

DIR=SGRT<DIRX*DIRX+DIRZ*DIRZ) 

IF(DIR*GE*1 *E-6)THEN 

UXM1=UX 

UYM1=UY 

UZM1=UZ 

UX=DIRX/DIR 

UY=0* ! FOR 2D CASE 5 UY*BIRY/DIR> FOR 3D CASE 

UZ=DIRZ/DIR 
END IF 
C 

GO TO 2000 
C 

3000 NFL=NL 

IF(K*LT*1*0R*K*GT»KT0T*QR» 

+ RX(NL-l) *LT*0*OR*RX(NL-1) ♦GT*XTQT)THEN 
WRI TE ( 2 * * ) ' LEADER OUT OF BOUNDS ' 

END IF 

IF ( EF AMB . LT ► BDF ) THEN 

URITE<2>») 'LEADER SEGMENT ENTERED REGION LT BREAK DOWN' 

END IF 
C 

CLOSE (1) 

CLOSE (2) 

RETURN 

END 

C 

C 

SUBROUTINE RANDOM (ISEEDt RAND) 

C RANDOM NUMBER GENERATOR 
C 

DOUBLE PRECISION E31*E32fSI* SJ 

E31»2***31 

E32=2.**32 

SI=DFLQAT ( ISEED ) 

SJ=MQD ( 69069 ♦fcSI+l * t E32) 

IF(ABS(SJ) *LT*E31 )THEN 
SI=SJ 

ELSEIF ( S J . GE ♦ E31 ) THEN 

SI=SJ-E32 

ELSE 

SI=SJ+E32 
END IF 
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IF(SI *GE»0* )THEN 

RAND=SI/E32 

ELSE 

RAND=1*+SI/E32 
END IF 

ISEED=INT(SI) 

RETURN 

END 



LDPM.FOR 
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C DEVELOPED BY E*M*A* - AUGUST 1986 

C PROGRAM DETERMINES DIRECTION OF NEXT LEADER SEGMENT DUE TO 

C ELECTRON WHICH WILL FORM AN IONIZED PATH TO THE SEGMENT TIP IN 

C THE LEAST TIME 

E FIELD TOTAL * E FIELD SEGMENT + E FIELD PREVIOUS SEGMENT 
+ AMBIENT E FIELD 

NCEL- # CELLS IN X AND IN Z 
ICXfICZ* XfY COORD OF INITIAL CHARGE REGION 
PARAMETER ( NCELX=97 f NCELZ=97 f ICX=8A f ICZ=31 ) 

DIMENSION FEX < 9409 )t FEZ (9409) 

CHARACTERS INFILE f0FILE1fGFILE3 
C 

C ISEED=RANDQM FUNCTION INTEGER 
ISEED=20897 

C SPECIFIED CHARGED REGION CELLS (OTHER THAN INITIAL CHARGED CELL) 
DATA NC1fNC2/195f18/ 

C 

INFILE='EFS*DT' ! INPUT AMBIENT E FIELDS 

QFILEi*'LDPM3*DAT' ! BINARY PLOTFILE OF LIGHTNING PATH 

QFILE3='LDPM.PR1NT' ! PRINT VALUES OF LIGHTNING MODEL 
C 

C INPUT VARIABLES 

C CL=SEGHENT LENGTH ELAMBDA= SEGMENT CHARGE DENSITY 

C BDF=BREAK DOWN FIELD EK=l/4PIEf>so NL=LEABER SEGMENT * 

C DZ=Z CELL HEIGHT DX=X CELL WIDTH 

C NE=# FREE ELECTRONS ED=ELECTRON DENSITY (.2/»«3) 

C Q=MAX ALLOWABLE ANGLE BTW EFIELD 2 DIRECTION VECTORS 
C 

DATA CL > ELAMBDA f BDF/50 * r * 001 f 1 . 5E5/ 

DATA PIfEKfNL/3*1415927f9*E9f0/ 

DATA DZfDXfNEfED»Q/200»f200»f500f »2f12»5/ 

DATA EXM1 FEYMi»EZMl/3*0*/ 

TWGPI=2>*PI 

QP=Q*PI/180. 

CONST =EK*EL AHBDA 
VOL=FLOAT (NE)/ED 
RADIUS** »75*VQL/PI>«<l,/3.) 

KTOT=NCELX*NCELZ 

XTOT=NCELX*DX 

C 

C READ IN AMBIENT E FIELDS 

OPEN ( 2 f FILE=INFILEf FORM* ' UNFORMATTED ' f STATUS* 'OLD' ) 
READ(2)(FEX<K)fK*1fKTOT)f(FEZ(K)fK*1fKTOT) 

CLOSE (2) 

C 

C INITIALIZE LEADER SEGMENT PROPAGATION 
C INPUT DATA: R1f AMBIENT E FIELDf U1 
C COORDINATES OF INITIAL POINT ON GRID PROBLEM SPACE 
RX*DX*(FLQAT<ICX)-.5) 

RY=0# 

RZ=DZ*<FLQAT(ICZ>-.5) 

C 

IC=ICX+NCELX*( ICZ-1 ) 

ETOT=SQRT ( FEX ( IC ) «2+FEZ ( IC ) **2 ) 

C DEFINE UNIT VECTOR 
UX=FEX< IO/ETOT 
UY=0* 

UZ=FEZ( IO/ETOT 
C 
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GPEN< 1 * F I LE=0FILE1 1 FORH* ' UNFORMATTED ' , ST ATUS* 7 NEW ' ) 

OPEN ( 3 , F ILE=QFILE3 r STATUS 3 ' NEW 9 ) 

C 

C START LOOP 
2000 NL*NL+1 

DETERMINE CELL POSITION OF NEW LEADER SEGMENT TIP ON OVERLAY 
USE APPROPRIATE AMBIENT E FIELD MAGNITUDE AND DIRECTION 
K IS THE CELL NUMBER OF THE LEADER SEGMENT TIP LOCATION 
KX=ANI NT ( ( RX/DX ) + * 5 ) 

KZ=ANINT ( <RZ/DZ)+*5) 

K=KX+NCELX*<KZ-1) 

EF AMB=SQRT < FEX ( K ) W2+FEZ ( K ) **2 ) 

TERMINATE IF 

LEADER SEGMENT HAS PROPAGATED OUTSIDE THE GRID 
LEADER SEGMENT ENTERING CELL WITH AN AMBIENT FIELD LT 1*5M V/M 
LEADER SEGMENT ENTERS INTO CHARGED REGION 
IF(K*LT*1 *QR*K*GT*KTOT*GR* 

+ RX*LT.0*0R*RX.GT*XTQT)6QT0 3000 

WRITE(3,*)NL,K,KX,KZiRX,RZ»ET0T?EFANB,FEX(K)»FEZ(K) 

PRINT*, NL,K,KX,KZfRX,RZf£TOT,EFAMB,FEX<K),FEZ(K) 
WRITE(1)RX,RZ 


IF(EFAMB*LT »BDF)GOTO 3000 

SHORTEST TIME PATH OF ELECTRON TO LEADER SEGMENT TIP 

THIS TIME PATH DETERMINES DIRECTION OF NEW LEADER SEGMENT 
TMIN=1*E20 

USE MONTE CARLO METHOD FOR FREE ELECTRONS 

TO CALCULATE ELECTRON DISTRIBUTION NEAR LEADER SEGMENT TIP 
RXP=RX+UX*CL 
RYP=RY+UY*CL 
RZP=RZ+UZ*CL 

DO 100 1=1, NE 

IN WORLD COORDINATES 

CALL RANDOM ( ISEED,RAND) 

RE=RAND*RADIUS 
RE=SQRT(RAND*(RADIUS**2) ) 

CALL RANDOM < I SEED , RAND ) 

TH=RAND*PI 
TH=ACOS ( RAND*2-1 ) 

CALL RANDOMUSEEDrRAND) 

PH=RAND*TWOPI 
RE-RAN ( I SEED) *RADIUS 
TH=RAN ( ISEED)*PI 
PH=RAN<ISEED)*2.*PI 
RMX=RE*SIN(TH)*COS(PH) 

RMY=RE*SIN<TH)*SIN(PH) 

RMZ=RE*COS(TH) 

X=RMX+RXP 

Y=RMY+RYP 

Z=RMZ+RZP 

RM=SQRT < RHX**2+RMY«2+RHZt*2 ) 
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C COMPUTE E FIELD COMPONENTS FOR Ith ELECTRON 
AX*X-RX 
AY=Y-RY 
AZ-Z-RZ 

ASSORT ( AX**2t A Y**2+AZ*<2 ) 

IF(A.LT.l.E-<S)GOTO 100 

B-AX*UX+AY*UY+AZ*UZ 

AS=A**2 

BS a B**2 

IF(AS-BS*LT*1 *E-6>60T0 100 
P a CQNST / ( AS-BS ) 

D a CL**2-2**B*CLFAS 
IF<D»LE»0» )GOTO 100 
D*SQRT(D) 

C 

C E FIELD COMPONENTS DUE TO LEADER SEGMENT AND AMBIENT FIELD 
EX a P*< ( (AX*(CL-B)+UX*<AS-B*CL> )/D)- 
+ ( ( UX*AS-AX*B ) / A ) ) +FEX ( K ) 

E Y=P* < ( ( AY* ( CL-B ) +UY* ( AS-B*CL > > /D > - 
+ ( (UY*AS-AY*B)/A) ) !+FEY(K) 

EZ-P*( ( ( AZ*(CL-B)+UZ*(AS-B*CL) >/D>- 
+ ( (UZ*AS-AZ*B)/A> )+FEZ(K) 

C 

IF(NL.LE*1)G0TQ 1000 
AXM1=X-RXM1 
AYM1=Y-RYM1 
AZM1=Z-RZM1 

AMI =SQRT ( AXM1**2+AYM1«2+AZM1**2) 

IF(AM1.LT*1*E-6>GQTQ 100 

BM1-AXM1*UXM1+AYM1*UYM1FAZM1*UZM1 

ASM1=AM1**2 

BSM1=BM1**2 

IF( ASM1-BSM1 ♦LT»1*E-6)G0TO 100 
PM1=C0NST/(ASM1-BSM1) 

DM1 =CL**2-2 * *BM1*CL+ASH1 
IF < DM1 >LE.O* )GOTQ 100 
Dttl=SQRT<DMl) 

C 

C E FIELD COMPONENTS DUE TO PREVIOUS LEADER SEGMENT 

EXM1-PM1*<((AXM1*<CL-BM1)+UXM1*<ASM1-BM1*CL))/DM1)- 
+ ( (UXM1*ASM1-AXM1*BM1)/AM1) ) 

EYM1=PM1*< < ( AYM1*(CL-BH1 HUYM1*< ASM1-BM1*CL> >/DMl )- 
+ ( ( U YM 1 *ASM 1 -AYMKBhl ) /AMI ) ) 

EZM1=PM1*(((AZM1*(CL-BH1)+UZM1*<ASM1-BM1*CL))/DMD- 
+ ( (UZM1*ASM1-AZM1*BM1)/AM1 ) ) 

C 

C TOTAL E FIELD COMPONENTS FOR Ith ELECTRON 
1000 EX=EX+EXMl 
EY=EY+EYM1 
EZ-EZFEZM1 

ET=SQRT ( EX**2+EY**2+EZ«2 ) + ♦ 01 
BL= ( RMX*EX+RMY*EY+RMZ*EZ ) / ( RM*ET ) 

C 

C CHOOSE APPROPRIATE ELECTRON TO BE DRAWN TOWARD 

C LEADER SEGMENT TIP IN SHORTEST TIME 

IF(ET *LT*(3*0E6) *AND* (RM/ET ) »LT* TMIN * 

F AND * ACOS < ABS ' DL ) ) ♦ LE * QP ) THEN 
TMIN=RM/ET 
ETOT=ET 
DIRX=EX 
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DIRY=EY 
DIRZ-E2 
END IF 

100 CONTINUE 
C 

C CALCULATE DIRECTION FOR NEXT LEADER SEGMENT 
DRX=CL*UX 
DRY=CL*UY 
DRZ=CL*UZ 
RXM1=RX 
RYN1*RY 
RZMZ=RZ 
RX=RXFDRX 
RY=RY+DRY 
RZ=RZ+DRZ 
C 

C NEW UNIT VECTOR 

C DIR=S0RT<DIRX**2+DIRY**2+DIRZ**2) 

DIR=SQRT<DIRX**2+DIRZ**2) ! PROJ'N ONTO 2D GRID 
IF(DIR*LT*1 *E-6)6QT0 2000 
UXM1=UX 
UYM1=UY 
UZM1=UZ 
UX=DIKX/DIR 
C UY=DIRY/DIR 

UY=Q. 

UZ=DIRZ/DIR 

C 

GO TO 2000 
C 

3000 IF(K*LT»1.QR*K*GT»KTOT»OR» 

+ RX.LT .0*0R*RX*GT»XTQT)THEN 

PRINT** 7 LEADER SEGMENT OUT OF BOUNDS' 

WRITE(3f*) 'LEADER SEGMENT OUT OF BOUNDS' 

END IF 

IF < EFAMB.LT .BDF)THEN 

PRINT* f 'LEADER SEGMENT ENTERED REGION LT BREAK DOWN' 
WRITE(3f*> 'LEADER SEGMENT ENTERED REGION LT BREAK DOWN' 
END IF 
C 

CLOSE ( 1 ) 

CLQSE<3> 

STOP 

END 

C 

C 

SUBROUTINE RANDQM( ISEEDrRAND) 

RANDOM NUMBER GENERATOR 

DOUBLE PRECISION E31>E32*SIrSJ 
£31=2. **31 
E32=2.**32 
SI=DFLOAT ( ISEED) 

S J=MOD < 6906? . *SI +1 . , E32 ) 

IF(ABS(SJ> *LT .E3DTHEN 
SI=SJ 

ELSEIF(SJ.GE*E31)THEN 
5I=SJ-E32 
ELSE 

SI=SJ+E32 
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END IF 

IF<SI *GE*0* )THEN 

RAND a SI/E32 

ELSE 

RAND-1 *+SI/E32 
END IF 

ISEED*INT(SI) 

RETURN 

END 



LDPM-INTERACTIVE 
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INTERACTIVE VERSION OF LDP H 


PROGRAM DETERMINES DIRECTION OF NEXT LEADER SEGMENT DUE TO 
ELECTRON WHICH WILL FORM AN IONIZED PATH TO 
THE STEPPED LEADER SEGMENT TIP IN THE LEAST TIME 
E FIELD TOTAL = E FIELD LEADER + E FIELD PREVIOUS LEADER 
+ AMBIENT E FIELD 


C 


C 

c 


c 


c 


10 

11 

c 

101 


CQMM0N/MAIN/FEX(?409) ? FEZ(9409) 

CQMM0N/VAR/NCX(10) ?NCZ(10) ?NCK<10) ?QN( 10) 
COMMON/DAT/NCELX?NCELZ?ICX?ICZ?NC?QP 
CHARACTER*40 INFILE ?QFILE1 ?0FILE3 

0FILE1 - ' LDPM.DAT ' ! BINARY PLOTFILE OF LIGHTNING PATH 

QFILE3='LDPN* PRINT' ! PRINT VALUES OF LIGHTNING MODEL 


PRINT*? 'LIGHTNING PROPAGATION - TORTUOSITY MODEL' 

PRINT*? 'DEVELOPED BY E.M.A. - AUGUST 1986' 

PRINT*? ' 

PRINT*? 'SPECIFY GRID INPUT DATA' 

PRINT*? ' 

PRINT*? 'GRID SIZE - ENTER THE TOTAL NUMBER OF CELLS' 
PRINT*? 'EACH CELL IS 200* >i 200*' 

PRINT*? 'ALONG THE X (HORZ AXIS)? ALONG THE Z (VERT AXIS)' 
PRINT*? 'ENTER X?Z (MAX IS 97?97)' 

READ*? NCELX?NCELZ 

PRINT*? ' 

PRINT*? 'ENTER THE INITIAL CELL COORD (X?Z) OF DISCHARGE' 
PRINT*? 'ENTER X?Z' 

READ*? ICX?ICZ 

PRINT*? 'ENTER THE CHARGE ( + ?-COUL ) ' 

READ*? QP 

PRINT*? ' 

PRINT*? 'YOU MAY SPECIFY ANY CHARGE REGION COORDINATE ' 
PRINT*? 'OTHER THAN THE INITIAL DISCHARGE COORDINATE. ' 
PRINT*? 'IF YOU DONT WANT THESE REGIONS ENTER 0 ' 

PRINT*? 'IT IS SUGGESTED THAT OVER THE ENTIRE GRID ' 
PRINT*? 'THE SUM OF THE TOTAL CHARGE = 0 ' 

PRINT*? 'ENTER THE # OF CHARGE REGIONS (MAX IS 10)' 
PRINT*? '(NOT INCLUDING THE INITIAL DISCHARGE COORDINATE >' 
READ*? NC 

IF(NC.EQ.0)G0T0 101 

PRINT*? 'ENTER THE X?Z PAIRS AND THEIR CHARGE ( - ? +COUL ) ' 
PRINT*? 'ie; XI ?Z1 ?CHARGE1 ' 

PRINT*? ' X2?Z2?CHARGE2' 

DO 10 I=1?NC 

READ*? NCX(I) ?NCZ( I ) ?QN(I) 

DO 11 1=1 ?NC 

NCK ( I ) =NCX ( I ) +NCELX* (NCZ(I)-l) 

PRINT*? ' 

PRINT*? 'THE PROGRAM USES A MONTE CARLO METHOD TO' 

PRINT*? 'GENERATE FREE ELECTRONS AROUND EACH STEPPED' 
PRINT*? 'LEADER TIP. AN ASSUMED DENSITY IS .2frlec/»**3' 
PRINT*? '300 ELECTRONS IS RECOMMENDED' 

PRINT*. 'LESS THAN 100 WILL LIMIT TORTUOSITY' 

PRINT*? 'ENTER THE NUMBER OF ELECTRONS DESIRED' 
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READ** NE 
C 

PRINT** ' 

PRINT** 'THE MONTE CARLO ROUTINE REQUIRES A SEED INTEGER' 
PRINT** 'ENTER ANY INTEGER UP TO 5 DIGITS l LT 32000.' 
READ** ISEED 
C 

PRINT** ' 

PRINT** 'PROGRAM GENERATES A PRINT FILE OF THE MODEL' 
PRINT** 'VALUES I A BINARY PLOTFILE IN AN X*Z OUTPUT' 
PRINT** 'AMBIENT E FIELDS ARE PRINTED OUT' 

PRINT** 'THESE FILES ARE LDPM .PRINT? LDPM.DAT* AMB.DAT' 
PRINT** ' 

PRINT* * ' *************************************** ' 

PRINT** ' 

C 

C 

C CALCULATE AMBIENT E FIELDS 
CALL EGRID 

C INPUT VARIABLES 

C CL=LEADER LENGTH ELAMBDA=LEADER CHARGE DENSITY 

C BDF=BREAK DOWN FIELD EK=1/4PIEpso NL=LEADER * 

C DZ=Z CELL HEIGHT DX=X CELL WIDTH 

C NE=* FREE ELECTRONS ED=ELECTRON DENSITY (.2/a**3) 

C Q=MAX ALLOWABLE ANGLE BTW EFIELD l DIRECTION VECTORS 
C 

DATA CL * EL AMBDA , BDF/50 .*,001*1. 5E5/ 

DATA PI *EK*NL/3, 1415927 r9.E9*0/ 

DATA BZ*DX*NE*ED*Q/200. *200. *500* .2*12.5/ 

DATA EXM1 *EYM1 *EZMl/3*0./ 

CONST =EK*EL AMBDA 
VOL=FLOAT(NE)/ED 
RADIUS=( .75*V0L/PI >**< 1 ./3. > 

KTQT-NCELX*NCELZ 

XTOT=NCELX*BX 

C 

C INITIALIZE LEADER PROPAGATION 
C INPUT DATA: Rlr AMBIENT E FIELDf U1 
C COORDINATES OF INITIAL POINT ON GRID PROBLEM SPACE 
RX=DX* ( FLOAT < ICX)- *5) 

RY=0» 

RZ=DZ*<FL0AT<ICZ>-.5> 

C 

IC=ICX+NCELX*( ICZ-1 ) 

ETOT =SQRT ( FEX ( I C ) **2+FEZ ( IC) **2 ) 

C DEFINE UNIT VECTOR 
UX=FEX< IO/ETOT 
UY=0. 

UZ=FEZ( IO/ETOT 
C 

OPEN ( 1 * FILE=OF I LEI f FORM 3 ' UNFORMATTED ' fSTATUS-'NEW' ) 

OPEN < 3 *FILE=0FILE3* STATUS* 'NEW' ) 

C 

C START LOOP 
2000 NL=NLF1 
C 

C DETERMINE CELL POSITION OF NEW LEADER SEGMENT TIP ON OVERLAY 

C USE APPROPRIATE AMBIENT E FIELD MAGNITUDE AND DIRECTION 

C K IS THE CELL NUMBER OF THE LEADER SEGMENT TIP LOCATION 
KX=ANINT ( (RX/DXH.5) 
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KZ-ANINT ( ( RZ/DZ ) + ♦ 5 ) 

K=KXfNCELX*(KZ-l> 

EFAMB=SQRT ( FEX < K ) **2+FEZ (K)**2) 

C 

C TERMINATE IF 

C LEADER HAS PROPAGATED OUTSIDE THE GRID 
C LEADER ENTERING CELL WITH AN AMBIENT FIELD LT I >511 V/H 
C LEADER ENTERS INTO CHARGED REGION 
IF(K.LT>l*OR*K*GT>KTQT*QR* 

+ RX*LT*0>QR,RX*GT?XTGT)6GT0 3000 
C 

WRITE(3?*) 'LEADER SEGHENT NUMBER '?NL 
WRITE<3?*)'CELL NUMBER '? K?' <X?Z> CQQRB'fKXfKZ 

WRITE(3?*) 'TOT E FIELD V/a'rETOTr' AMB E FIELD ' fEFAHB 
URITE(3? *) 'AMB EX 'fFEX(K)?' ANB EZ 'fFEZ(K) 

WRITE<3?*) 'X *et> frcm 0?0 '?RXf' Z *et* tram OfO 'fRZ 
WRITE(3?*>' 

PRINT*? 'LEADER SEGMENT NUMBER '?NL 
PRINT*? 'CELL NUMBER 'fK?' <X?Z> COORD' fKXfKZ 

PRINT*? 'TOT E FIELD V/«'?ETOT?' AMB E FIELD 'fEFAHB 
PRINT*? 'AMB EX '?FEX(K)?' AMB EZ 'tFEZ(K) 

PRINT*? 'X met* fr on 0?0'?RX?' Z »et* from 0?0 '?RZ 
PRINT*? ' 

URITE<1 )RX?RZ 
C 

IF(EFAMB.LT *B0F)GQTQ 3000 
IF(NC*EQ*0* >80T0 55 
DO HI II-1?NC 

111 IF(K>£Q*NCK(II) )GOTQ 3000 
C 

C SHORTEST TIME OF ELECTRON TO FORM IONIZATION CHANNEL 
C THIS TIME DETERMINES DIRECTION OF NEW LEADER SEGMENT 
55 TMIN«1 >E20 

C 

C USE MONTE CARLO METHOD TO GENERATE RANDOM FREE ELECTRONS 
C CALCULATE ELECTRON POSITIONS NEAR LEADER TIP 
DO 100 I-1?NE 

IN WORLD COORDINATES 

CALL RANDOM < I SEED? RAND) 

RE"RANO*RADIUS 
C RE-SQRT < RAND* ( RAD IUS**2 ) ) 

CALL RANDOM ( ISEED ? RAND ) 

TH=(RAND)*PI 

C TH=ACQS<(RAND*2>-1) 

CALL RANDOM (ISEED? RAND) 

PH=(RAND)*2.*PI 

X«RE*S IN ( TH ) *COS < PH ) +RX+UX*CL 

Y=RE*S I N ( TH ) *S I N ( PH ) +R Y+UY*CL 

Z=R£*CQS ( TH ) +RZ+UZ*CL 

RMX=X-RX-UX*CL 

RMY*Y-RY-UY*CL 

RHZ S Z-RZ-UZ*CL 

RM=SQRT ( RMX**2FRMY**2+RMZ**2 ) 

C 

C COMPUTE E FIELD COMPONENTS FOR Ith ELECTRON 
AX=X-RX 
AY=Y-RY 
AZ=Z-RZ 

ASSORT ( AX**2+A Y**2+AZ**2 ) 
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IF( A*LT*1 *E-6)GQTQ 100 

B=AX*UX+AY*UY+AZ*UZ 

AS=A**2 

BS=B**2 

IF(AS-8S.LT.1*E-6)G0TQ 100 
P=CGNST/<AS-BS) 

D=CL**2-2.*B*CL+AS 
IF(D*LE*0* )GQTQ 100 
D=SQRT(D) 

C 

C E FIELD COMPONENTS DUE TO LEADER SEGMENT AND AMBIENT FIELD 
EX=P*( ( ( AX*(CL-B)+UX*<AS-B*CL> )/D>- 
+ ( ( UX* AS- AX*B )/A) ) +FEX ( K ) 

EY=P*< < <AY*(CL-B)+UY*(AS-B*CL) >/D>- 
+ (<UY*AS-AY*B)/A)) !+FEY(K) 

EZ=P* < ( ( AZ* ( CL-B ) HIZ* ( AS-B*CL ) ) /D ) - 
+ ( ( UZtAS- AZ*B )/A) ) +FE2 ( K ) 

C 

IF(NL*LE»1)G0T0 1000 
AXN1»X-RXM1 
AYM1=Y-RYM1 
AZM1=Z-RZM1 

AMl=SQRT(AXMm2+AYMl*mAZHl«2) 

IF (AMI ♦LT»1*E-£)GQTG 100 

BM1=AXM1*UXM1+AYM1*UYM1+AZM1*UZM1 

ASM1=AM1**2 

BSM1=BM1**2 

IF(ASM1-BSM1 .LT*1*E-6)G0T0 100 
PMl=CGNST/< ASH1-BSM1 ) 

DM1=CL«2-2>*BM1*CL+ASH1 
IF (DM1 #LE*0*)GQTQ 100 
DM1»S0RT(DM1) 

C 

C E FIELD COMPONENTS DUE TO PREVIOUS LEADER SEGMENT 

EXM1=PM1*(<<AXM1*<CL-BM1)+UXMU<ASM1-BM1*CL))/DM1>- 
+ ( (UXMlfASMl-AXMKBMl )/AMl ) ) 

EYM1*PM1*<(<AYN1*<CL-BH1)+UYN1*<ASH1-BM1*CL)>/DM1>- 
+ ( (UYMltASMl-AYMl<BMl)/AMl ) ) 

EZM1=PM1*< ( (AZMlt(CL-BMl )+UZMlt(ASMl-BMltCL) >/DMl)- 
+ ( (UZMltASMl-AZMl*BMl)/AMl) ) 

C 

C TOTAL E FIELD COMPONENTS FOR Ith ELECTRON 
1000 EX=EXFEXM1 
EY=EY+EYM1 
EZ=E2FE2M1 

ET=SQRT ( EX**2+EY**2+EZ«2 ) + . 01 
DL= ( RHX*EX+RMY*EY+RMZ*EZ ) / ( RH*ET ) 

C 

C CHOOSE APPROPRIATE ELECTRON TO FORM IONIZED CHANNEL 
C TO THE LEADER TIP IN THE SHORTEST TIME 

IF(ET*L?.(3.0E6KAND.<RH/ETKLT*TMIM. 

+ AND*AGGS(ABS(DL) ) *LE.Q*PI/180* )THEN 
TMIN=RM/ET 
ETQT=ET 
DIRX=EX 
DIRY=EY 
DIRZ=EZ 
END IF 

100 CONTINUE 
C 
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C CALCULATE DIRECTION FOR NEXT LEADER SEGMENT 
DRX=CUUX 
DRY=CL*UY 
DRZ=CL*UZ 
RXM1=RX 
RYM1-RY 
RZMZ=RZ 
RX-RX+DRX 
RY=RY+DRY 
RZ=RZ+DRZ 

NEU UNIT VECTOR 

DIR=SQRT ( DIRXW2+DIRY**2+DIRZ**2 ) 
DIR«SQRT<DIRX**2+DIRZ**2> ! PRQJ'N ONTO 2D GRID 
IF<DIR*LT.1,E-6)GOTO 2000 
UXM1 S UX 
UYM1 3 UY 
UZM1=UZ 
UX=DIRX/DIR 
C UY=DIRY/DIR 

UY=0* 

UZ=DIRZ/DIR 

C 

GO TO 2000 
C 

3000 IF<K*LT*1*0R*K*GT *KTGT*OR* 

+ RX *LT*0*0R*RX*GT ♦XTQTJTHEN 
PRINT*? 'LEADER OUT OF BOUNDS' 

WRITE < 3?*) 'LEADER OUT OF BOUNDS' 

END IF 

IF<EFAHB.LT*BDF)THEN 

PRINT*? 'LEADER ENTERED REGION LT BREAK DOUN' 
URITE< 3?* >' LEADER ENTERED REGION LT BREAK DOUN' 
END IF 

DO 112 II a l>NC 

IF<K*EG*NCK(II) )THEN 

PRINT*? 'LEADER ENTERED CHARGED REGION' 

WRITE<3?*> 'LEADER ENTERED CHARGED REGION' 

END IF 

112 CONTINUE 
C 

CLOSE (1) 

CL0SE(3) 

STOP 

END 

C 

C 

SUBROUTINE RANDOM ( I SEED? RAND) 

C RANDOM NUMBER GENERATOR 
C 

DOUBLE PRECISION E31 ?E32»SI?SJ 

E31=2***31 

E32=2.**32 

SI=DFLOAT( ISEED) 

SJ=MOD< 69069 ♦*SI+1 » ?E32) 

IF(ABS(SJ) *LT *E31)THEN 
SI=SJ 

ELSEIF(SJ*GE»E31)THEN 

SI=SJ-E32 

ELSE 
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SI=SJ+E32 
END IF 

IF(SI »GE»0»)THEN 

RAND=SI/E32 

ELSE 

RAND=1 *+SI/E32 
END IF 

ISEED=INT(SI) 

RETURN 

END 

C 

C 

SUBROUTINE E6RID 

C CALCULATION OF AMBIENT ELECTRIC FIELDS 
C 

DIMENSION X<9409>»Z<9409> 

DIMENSION XMRN<5)fB<5> 

COMMON/MAIN/FEX ( 9409 ) , FEZ ( 9409 ) 
COMMON/VAR/NCX(10)fNCZ(10)fNCK(10)fQN(10) 

COMMON/DAT /NCELX f NCELZ > ICX » ICZ r NC t QP 
DATA DXrDZ*DP>PI/200t >200 » »9*E9 f 3* 1413927/ 

C 

OPEN ( 4 * FI LE= ' AMB * DAT ' » ST ATUS* ' NEW ' ) 

C 

DO 10 >1, NCELZ 
DO 10 1=1 , NCELX 
K=I+NCELX*(J-1) 

X<K>=DX*<FLQAT<I>-*5> 

Z(K)=DZ*CFL0AT(J)-,3) 

K1=K 

10 CONTINUE 
C 

LP=ICX+NCELX*< ICZ-1 > 

C CALCULATE RESULTANT E FIELDS FROM CHARGE CENTERS 
DO 20 K*i>Kl 
IF<NC*EQ*0>THEN 

IF(X(K) *EQ.X(LP) »AND*Z(K) *EQ*Z(LP) )GOTO 20 
END IF 

IF(NC»EQ»0)G0T0 23 
DO 333 1=1 fNC 

IF(X(K) *EQ»X(LP) »AND»ZOO *EO.Z(LP) *0R* 

+ X(K) »EQ»X(NCK<I) )*AND*Z(K)»EQ»Z(NCK(I)))GOTO 20 
333 CONTINUE 

C 

23 XMRP=SQRT<(X(K)-X(LP))«2+(Z(K>-Z(LP)>«2) 
A=DP*QP/XMRP**3 

FEX(K)=0* 

FEZ(K)=0> 

IF<NC*EQ*0)G0TQ 24 
DO 444 1=1 »NC 

XMRN(I)=SQRT((X(K)-X(NCK(I)))«2+(Z(K)-Z(NCK(I)>>«2> 
B < I ) =DP*QN ( I ) /XMRN ( I ) «3 
FEXCK)=FEX<K>+BCI)*(X<K)-X<NCK<I>)> 
FEZ(K)=FEZ<KHB(I)*<Z<K>-Z(NCK<I>>> 

444 CONTINUE 

C 

24 FEX(K)=FEX<KHA«X<K)-X(LP>> 
FEZ(K)=FEZ(K)+A*<Z(K)-Z(LP> ) 

20 CONTINUE 

FEX ( LP ) -FEX ( LP+1 ) 
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FEZ ( LP ) =FEZ ( LPM > 

IF(NC»EQ»0)6QT0 25 
DO 555 1*1 iNC 
FEX(NCK(I) )=FEX(NCK< I)-l) 

FEZ(NCK(I) )*FEZ(MCK(I)-1) 

555 CONTINUE 

25 DO 26 K=1»K1 
WRITE(4p<)K>FEX(K) fFEZ(K) fX(K) fZ(K> 

26 CONTINUE 
CLOSE < 4) 

RETURN 

END 
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